Molecular Binding Energies from Partition Density Functional Theory 
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Approximate molecular calculations via standard Kohn-Sham Density Functional Theory are 
exactly reproduced by performing self-consistent calculations on isolated fragments via Partition 
Density Functional Theory [Phys. Rev. A 82, 024501 (2010)]. We illustrate this with the binding 
curves of small diatomic molecules. We find that partition energies are in all cases qualitatively 
similar and numerically close to actual binding energies. We discuss qualitative features of the 
associated partition potentials. 



I. INTRODUCTION 

Kohn-Sham Density Functional Theory (KS-DFT) 
[TJ [5] provides one of the most useful methods for cal- 
culating electronic properties of molecules and materials 
[3]. The accuracy that can be achieved with modern ap- 
proximations to the xc- functional [IHS] , and the efficiency 
of numerical implemenations [7], make of KS-DFT one 
of the workhorses of computational quantum chemistry, 
materials science, and nanotechnology. 

A method has recently been proposed that exactly 
reproduces the results of approximate KS-DFT calcu- 
lations via self-consistent calculations on isolated frag- 
ments '8]. Partition Density Functional Theory (PDFT) 
is based on the density- Partition Theory (PT) of ref . j9j , 
and is analogous to KS-DFT in that it establishes a map 
between the physical system of ground-state density ri(r), 
thought of as a collection of interacting fragments, and 
an auxiliary system of Nf non-interacting fragments with 
ensemble-ground-state-densities {^^(r)} subject to the 
density constraint X)q^ '^a(r) = n{Y). The appeal of this 
formalism is two-fold: on the one hand, by optimally di- 
viding a complex system into fragments it allows one to 
build a rigorous foundation for chemical reactivity the- 
ory [SHU]- On the other hand, by solving the molecular 
Kohn-Sham equations in a different way, it focuses at- 
tention on quantities that are amenable to new, differ- 
ent approximations, potentially leading to linear-scaling 
algorithms for large systems. In that spirit, PDFT is 
similar to embedded-DFT [T^15j . One such quantity is 
the partition potential, Vp{r), a global molecular prop- 
erty (called reactivity potential in ref. ^9J, and equivalent 
in practice to the crystal potential of ref. [IS] and the 
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embedding potential of ref.[T3]). It is used as an addi- 
tional external potential for each fragment so that the 
density constraint is satisfied. If the density of the whole 
is known, as in the original formulation of PT, Up(r) is 
simply the Langrage multiplier of the density constraint. 
Without knowing the whole density, an iterative pro- 
cedure for deriving Vp{r) was designed and illustrated 
with 1-D models [HI [T7] . In this work we will show how 
this procedure is realized for real molecules. This step is 
of course essential to be able to explore the promise of 
PDFT on both fronts mentioned before. 

The partition potential is the functional derivative of 
the partition energy, Ep[{na}], with respect to any of 
the fragment densities na(r). This energy is defined as 
the difference between the molecular ground-state energy 
E[n\ and the sum of the fragment energies X^a^ Ea[na\. 

Kohn-Sham equations can be established so that ap- 
proximations to £'p[{nQ,}] yield definite predictions for 
the ground-state energy and density of the assembly. 
When the exact i?p[{7iQ}] is employed implicitly by in- 
version of molecular KS equations ^ , then the exact KS- 
DFT results are recovered without ever having to solve 
the direct problem for the assembly, but only for the frag- 
ments. 

Analytical studies on one-dimensional models of het- 
eronuclear diatomics [18J provide an early indication that 
the fragment dipoles obtained by PDFT are more ad- 
justed to chemical intuition and more transferable than 
those obtained by other density-partitioning schemes, 
but more studies are of course needed in real systems. 

In this work, by employing the Wu-Yang algorithm 
[19) for iterative inversion, we demonstrate convergence 
of the PDFT equations in small diatomic molecules, and 
discuss qualitative features of partition potentials and 
partition-energy binding curves for IIe2, II2, and LiH. 
We show that the partition energies and potentials are 
interesting quantities in themselves, as they can be used 
as conceptual and interpretative tools. 



First, we summarize the PDFT procedure in SecjTTj 
providing details of our implementation. Convergence of 
the PDFT equations is demonstrated in Sec |III| for the 
binding curves of He2, H2, and LiH, along with impli- 
cations, qualitative features of patition potentials, and 
-Ep-binding curves (in addition to actual binding curves) . 
Concluding remarks are given in Sec|IV[ 



II. METHOD 

For the simplicity of discussion, we consider a com- 
pound with only two parts {A and B), but the method is 
equally applicable to any number of fragments. We also 
limit ourselves to fragments with fixed integer number of 
electrons, as in related recent work on embedding-DFT 
[T^1 - [T5] , only briefly discussing the issue of chemical po- 
tential equalization and fractional electron numbers. 

In PDFT, the total energy is expressed as 



E[n] = EA[nA] + EBins] + Ep[nA, ub] 



(1) 



where n(r) — ^^^(r) + ns(r), and a common functional 
for E, Ea and Eb is assumed. The above equation can 
be viewed as a formal and exact definition of Ep. To 
minimize E by variations of fragments' densities, which 
are built from their own sets of orbitals, we have the 
following Kohn-Sham equations: 



V^ + Va{r) + Vp{r) + WHxc['^Q](r) 



^-" 'C(r)=erC(r) 

(2) 

Here, a is a fragment index, i.e A or B in this work. The 
partition potential Vp{r) is common to both fragments, 
thus has no a-index. 

Vp{r) could be derived explicitly if we knew the func- 
tional form of -E'p[{nQ}]. Without an expression for Ep 
as an explicit functional of the {«„}, it is also possible to 
derive Vp{r) through an iterative procedure, which was 
first proposed in ref. |8] and we reiterate here. 

Suppose that we are at the beginning of the k-th iter- 
ation. We obtain all fragment densities ria by solving 
Eq. [2j We then construct a total pro-molecule density 
as n('')(r) = ^n^ (r). Because the effect of Wp(r) is to 
make n(r) the same as the true ground-state density of 
the whole system ns(r), the difference between n*^'^)(r) 

and ns{r) should be used a guidance to update Vp (r). 
For that, we do a constrained search to find the energy 
of h'^'^K i.e. 



£:[n('=)] = min E[n]. 



(3) 



We employ the direct optimization algorithm of Wu 
and Yang [19 , as used in calculating the frozen density 
energy in a recently-developed density based energy de- 
composition analysis |31j . Thus we rewrite the above 
equation as 



^[nC^)] =E^[fi^''^]+E, 



HXCl"- 



- min {r,[*]+Ex[*]} 
(4) 



for a general hybrid functional, where E'xi^] represents 
a fraction of the HF exchange energy calcualted from a 
Slater determinant \1/ that is constrained to yield n'-''\ 
At the end of this minimization, the effective potential 
for the molecular Kohn-Sham orbitals is 

(5) 
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i;off(r) = Va{r) + Whxc[^^'''](i") - WAlr) , 

!\ {r\ is inst ttip T-ficrrfl.nP'P mnltinlipr rorrpsT 



where vx(j) is just the Lagrange multiplier corresponding 
to the density constraint and is expanded by a linear com- 
bination of atom-centered Gaussian functions. Because 
v\{r) is used to force the density of the whole system 
to be n^''\ its reverse should have the effect of making 
h{r) more like ns{r). That is: we can set Vp{r) — —v\{r) 
and start the next iteration of fragment calculations. In 
practice, we update Vp{r) as follows: 
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(6) 



where i is the iteration number, and is a damping fac- 
tor between and 1 used to control convergence. In our 
calculation, we have used 9 = 1 01 9 — 0.25. The con- 
vergency criterion we use is |i?[rifc] — -E[ns]| < e, where 
e = 10~^; this guarantees the converged energy is the 
same at the ground-state energy. The alternative choice 
of \E[nk]— E[nk-i] \ < £ gives essentially the same results. 



III. RESULTS 

We demonstrate our calculations of Vp(r) with three 
simple examples of diatomic systems: IIe2, II2, and LiH. 
In all calculations, Dunning's aug-cc-pvTz basis set is 
used for molecular orbitals. 

The counter-poise (CP) method is used to account 
for any Basis Set Superposition Error (BSSE). This ap- 
proach is crucial in PDFT since Vp (r) adds features to the 
fragment's effective potential directly at the location of 
the other atom, precisely where the ghost basis functions 
are added [32> 

The partition potential is expanded by atom-centered 
Gaussian functions, and each center has five s-type 
functions, with even-tempered exponents of 2",ri — 
0, ±2, ±4. In the following discussion, we will use sev- 
eral energy terms. Suppose E^ and E'^ are the en- 
ergies of the fragments with no influence of the par- 
tition potential; E^ and Eg are their energies with 
the converged partition potential; and Eab is the en- 
ergy of the compound. Therefore the binding energy is 
-Ebind = Eab — {Ea + Eg), and the partition energy is 
Ep — Eab — {Ea + E^). We also define the prepara- 
tion energy as Ep^^p ^ {E\ + E%)- {E°^ + E%), which 
is the energy increase associated with the deformation of 
fragments. Clearly, i^bind = £^prcp + Ep. 

A. Helium Dimer 

Rare-gas dimers are known to be weakly bound due to 
van der Waals interactions, which are not accurately cap- 



tured by most density-functional approximations. How- 
ever, because our procedure is general and independent 
of the exchange-correlation functional, it is not criti- 
cal to have the correct binding curve. Instead, for a 
clear demonstration, we use Hartree-Fock exchange only, 
which is known to be purely repulsive between nonpolar 
closed-shell systems. As shown in Fig. [TJ the binding en- 
ergy for He2 is all positive and increases rapidly when the 
internuclear distance is shortened. It also shows that the 
preparation energy is very small, which means the defor- 
mation in He atoms is small, as expected in this system, 
though it starts to grow when the atoms are too close to 
each other. The repulsive nature of the interaction means 
that electron densities are pushed away from each other 
when the two He atoms are in close contact. Thus the in- 
ternuclear region has a density decrease, as shown in Fig. 
[2J In PDFT, this density difference is achieved through 
deformation of each atom, due to the action of the par- 
tition potential. In Fig. [31 wc plot Vp along the internu- 
clear axis at a few representative internuclear distances. 
Clearly, Vp is most positive in the internuclear region, cor- 
responding to the density deficiency. The magnitude of 
Vp decreases as the internuclear distance increases, until 
to a point that no u„ is needed. 
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FIG. 2: The density differences in He2 as compared to the 
original atoms along the line through both nuclei. The total 
difference (dashed line) is the sum of the deformation in each 
atom (solid line). The nuclei coordinates are R = ±0.8 A. 
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FIG. 3: The partition potentials for He2 at different internu- 
clear distances. The nuclei are at ±_R. 



B. Hydrogen Molecule 



FIG. 1: The Hartree-Fock energies for He2 at different inter- 
nuclear distances. 



It is notable that there are significantly more oscilla- 
tions in the partition potential than in the density differ- 
ences. Some of the oscillations are physical. But there 
are at least two other possible reasons contributing to 
the oscillations in Vp. One is pathological with gaussian 
densities, as nicely explained by Schipper, Gritsenko and 
Baerends I33j. The other is numerical and due to the 
fact that we expand Vp in a finite basis set [33] . We have 
used a small number of functions so as to limit the oscil- 
lations caused by the expansion. However, we are unable 
to use non-gaussian densities yet, which makes it difficult 
to determine the nature of the oscillations. 



For the covalently bonded molecule H2, the natural 
choice of partition is to use two open-shell H (OSH) 
atoms. Because their spins are paired up in the molecule, 
we only consider the total charge density. Mathemat- 
ically one could also use half-occupied closed-shell H 
atoms (CSH) as the fragments, thus without polariz- 
ing the spin. We study the energetics of both parti- 
tions as a function of the internuclear distance, using the 
B3LYP approximation to the exchange-correlation func- 
tional. For the H2 molecule, we only consider restricted 
Kohn-Sham (RKS) calculations. It is well-known that a 
restricted calculation does poorly for large internuclear 
distances. The erroneous behavior is evident from the 
binding energy curve when the OSH atoms are used as 
the reference. As shown in Fig. [4| E^bind approaches a 
positive value instead of zero. On the other hand, when 
the CSH atoms are used as the reference, i?bind does go 
to zero. However, it becomes too large at the optimal 
bond length. The two binding curves are simply differ- 



ent by a constant shift, and this shift comes from the fact 
that OSH and CSH have different energies in the B3LYP 
approximation, while they should be degenerate with the 
exact functional 1351. 
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FIG. 4: Top: The B3LYP energies for H2 at different inter- 
atomic distances. The optimized bond length is _D = 0.743 
A. Bottom: The partition potentials for H2 at different inter- 
nuclear distances. 



In PDFT, the differences in the choice of fragments 
will not matter if the partition energy can compensate 
for the difference and yield identical total energy. In our 
case here, the two Ep are indeed quite different. How- 
ever, the two Ep curves differ more than by a simple 
constant shift. The non-uniform difference can be appre- 
ciated by comparing the preparation energies. i?prop of 
OSH fragments is smaller at short internuclear distances 
than that of CSH fragments. However, the latter goes to 
zero at long distances while the former does not. At long 
distances, a restricted H2 is essentially two half-occupied 
closed-shell H atoms, so the asymptotic behavior is not 
surprising. But it is interesting to see that at short dis- 
tances, the OSH fragments pay less penalty to make their 
densities resemble that of the molecule. 



C. Lithium Hydride 

As another example, we consider the heteronuclear 
LiH. Within the formal partition theory, there is a unique 
choice of the fragments, with their chemical potentials 



equilibrated. Achieving equilibration requires treating 
fragments with fractional number of electrons in the 
spirit of PPLB [3B,_. In that case, the number of elec- 
trons in a fragment is also a variable to be optimized. 
Because the partition potential will be different when 
the fragments change, the optimization of both the par- 
tition potential and the number of electrons is mutually 
dependent and has to be achieved simultaneously. We 
will treat this complexity in the future. In this work, we 
simply use fixed fragments and derive the corresponding 
partition potential. 

Without the optimal fragments, we consider all possi- 
ble partitions. For LiH, there are two possibilities. First, 
we use neutral atoms. Second, we use Li"*" and H~. We 
do the partition at the optimized internuclear distance of 
1.59073 A. For the neutral partition, Ep^^p = 0.053257 
a.u. and Ep = —0.146407 a.u. For the ionic partition, 
^prop = 0.034395 a.u. and Ep = -0.300161 a.u. The 
larger partition energy in the ionic case could be the re- 
sult of Coulomb attraction. However, the preparation 
energy is smaller for the ionic partition, suggesting the 
LiH bond is closer to an ionic bond than a covalent bond. 
What is surprising is that the partition potential for the 
ionic case looks much stronger than the neutral one (Fig. 
[5]), despite causing less distortion in fragment's energies. 
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FIG. 5: The B3LYP partition potentials for LiH. The Li atom 
is at a: = and H at 2; = 1.59073 A. 




IV. CONCLUDING REMARKS 

Without having to solve directly the KS equations 
for the total external potential, we have shown how the 
PDFT algorithm of ref. [8^ provides in practice the same 
answers via fragment-KS equations. In addition, this 
method yields fragment densities, fragment energies, and 
a partition potential that is shared by all fragments such 
that the sum of their densities reproduces the correct to- 
tal density. 

Although no physical meaning can be attached to a 
parition potential beyond the one implied by its defini- 
tion (i.e. that it is the potential common to all fragments 



such that the sum of the fragment densities equals the 
total molecular density) , some generic features of parti- 
tion potentials seem to go in line with chemical intuition: 
they are positive when the interaction between fragments 
is repulsive (case of He2 within Hartree-Fock), and their 
average magnitude is larger when the interaction between 
fragments is stronger. Similarly, the strength of the in- 
teraction between fragments is loosely measured by the 
magnitude of the partition energy. No such conclusion 
can be drawn for the preparation energy, however, as 
shown for the case of LiH where a somewhat larger prepa- 
ration energy is associated with a much smaller partition 
potential (neutral vs. ionic partition). But the prepara- 
tion energy can tell us about the character of the bond, 
an aspect that we plan to study further in future work. 
The case of LiH also highlights the need to go beyond 
integer numbers of electrons in our implementation of 
PDFT. 

PDFT calculations also allow us to look at the dis- 
sociation problem from a different angle. For example, 
we found that open-shell fragments in H2 are preferred 
at short inter-nuclear separations in the sense that they 
pay less penalty to make their densities resemble that 
of the molecule, but close-shell fragments are preferred 
at long separations. The respective preparation energies 
cross near the Coulson-Fischer point. 

Finally, we point out that from weak (He2) to rela- 



tively strong (H2) chemical bonds, partition energies are 
qualitatively similar to actual binding energies, and nu- 
merically close to them (i.e. preparation energies are 
small in the cases studied) . This similarity of Sp-curves 
to their corresponding binding curves suggests that ap- 
proximations of ii^p[{nQ,}] as explicit functionals of the 
{uai} might be very useful for practical computations. 
Not only would they provide a direct way to obtain the 
partition potentials by functional differentiation, circum- 
venting the need of expensive inversion steps; sensible 
approximations would also lead to energies that are close 
to actual binding energies. This is analogous to what 
happens in KS-DFT, whose success is largely due to the 
fact that the sum of KS orbital energies is typically close 
to actual ground-state energies in chemical applications. 
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